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Abstract 

The  main  focus  of  this  project  has  been  the  investigation  of  various  questions 
related  to  the  p—  and  h  —  p  versions  of  the  finite  element  method,  including  mixed 
methods.  These  new  versions  differ  form  the  classical  /i-version  where  the  degree  p 
of  polynomials  used  is  kept  fixed  (usually  p  =  1  or  2)  and  accuracy  is  achieved  by 
decreasing  the  mesh-size  h.  In  the  p-version,  a  fixed  mesh  with  constant  h  is  used 
and  p  is  increased  for  accuracy.  The  h-p  version  combines  the  two  approaches. 

We  have  studied  various  theoretical  and  computational  questions  related  to  these 
methods.  The  specific  areas  we  have  covered  are  as  follows: 

1.  Optimal  approximation  results  for  the  p-  and  h  —  p  versions. 

2.  Approximation  of  boundary  conditions  by  the  p  and  h  —  p  versions. 

3.  The  p  and  h  —  p  versions  for  mixed  methods. 

4.  The  p  and  h  —  p  versions  for  integral  equations. 

5.  An  implementation  of  the  p  and  h  —  p  versions  for  MODULEF. 

6.  Numerical  methods  for  a  class  of  reaction-diffusion  problems. 

Some  of  our  theoretical  results  have  already  been  incorporated  into  the  commercial 
h-p  version  code  PROBE  that  is  being  used  by  several  industries  in  the  U.S., 
including  those  involved  in  aeronautical  engineering. 
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1.  Introduction 


The  main  focus  of  this  project  has  been  the  investigation  of  various  questions  related 
to  the  p-  and  h-p  versions  of  the  finite  element  method  (FEM),  including  mixed  finite 
element  methods. 

The  p-  and  h-p  versions  are  relatively  new  methods,  the  first  theoretical  papers 
for  which  appeared  in  1981  [A4j,  [A3].  One  of  the  main  differences  between  these 
methods  and  the  standard  h—  version  (which  is  the  classical  form  of  the  FEM)  is  the 
use  of  higher  order  elements.  In  the  /i-version,  accuracy  is  achieved  by  decreasing  the 
mesh  size  h  while  keeping  the  degree  p  of  elements  used  fixed  at  a  low  level,  p  =  1,2. 
In  the  p-version,  the  mesh  is  kept  fixed  (i.e.  h  is  constant)  and  accuracy  is  achieved 
by  increasing  the  degree  p  of  elements.  The  h-p  version  uses  simultaneous  increase  of 
the  degrees  and  mesh  refinement  to  obtain  the  desired  accuracy. 

Early  theoretical  results  showed  that  for  smooth  solutions,  the  rates  of  convergence 
obtained  by  these  new  methods  were  not  bounded  by  the  degree  of  polynomials  used 
(as  in  the  /i-version)  and  consequently  could  be  significantly  higher.  Moreover,  the 
rates  of  convergence  obtained  for  solutions  of  linear  elliptic  problems  (where  ra  type 
singularities  arise  at  the  corners  of  the  domain)  were  also  established  to  be  higher 
(roughly  twice  that  obtained  from  the  /i-version).  However,  none  of  these  results 
established  the  optimal  rates  of  convergence  that  could  be  expected.  (Our  first  goal 
was  to  establish  such  rates). 

The  p-  and  h-p  versions  formed  the  basis  of  the  industrial  program  PROBE  [A10], 
developed  by  Noetic  Tech,  St.  Louis  [with  a  first  release  in  1985  and  a  second  one  in 
1986|.  This  program  implements  these  versions  for  two-dimensional  problems  (plane 
elasticity  and  heat  transfer).  In  a  relatively  short  period  of  time,  PROBE  has  become 
popular  in  several  industries  in  the  U.S.,  including  those  involved  in  aeronautical 
engineering.  One  reason  for  this  is  the  fact  that  the  p-  and  h-p  versions  are  the  most 
efficient  known  finite  element  methods  for  the  analysis  of  linear  elliptic  problems  over 
polygonal  domains  (see  [A2|  for  details  and  d<-s»  lipl'em  •>(  industrial  usage  feedback). 
Some  of  our  theoretical  results,  particularly  those  on  the  treatment  of  inhomogeneous 
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essential  boundary  conditions  by  the  p-version  have  found  direct  application  in  the 
implementation  of  I’ROBE. 

The  p-version  has  also  been  used  in  other  programs  like  FIESTA,  developed  at 
ISMES,  Bergamo,  Italy.  A  three-dimensional  p-version  code  called  STRIPE  [Al]  has 
been  developed  at  the  FFA  (Aeronautical  Research  Institute  of  Sweden)  and  can  cur¬ 
rently  be  implemented  on  CRAY  computers.  Moreover,  a  p-version  implementation 
for  MODULEF  [B16]  (INRIA,  France)  has  been  developed  under  this  project  and 
described  in  the  next  section. 

2.  Results  achieved 

In  this  section,  we  describe  the  results  achieved  under  the  grant  over  the  past  three 
years.  A  list  of  categories  under  which  results  have  been  obtained  may  be  found  on 
the  cover  page. 

2.1.  Optimal  Approximation  Results  for  the  p-  and  h  —  p  versions 

One  of  the  basic  pieces  of  desirable  information  for  any  finite  element  method  is  the 
asymptotic  rate  of  convergence  that  may  be  expected  and  whether  or  not  this  rate  is 
optimal.  For  the  standard  conforming  finite  element  method  based  on  the  p-version, 
previous  results  established  in  [A4]  by  Babuska,  Szabo  and  Katz  and  in  [A7|  by  Dorr 
were  non-optimal.  More  precisely,  for  a  two-dimensional  polygonal  region  H,  the  error 
e  in  approximating  a  function  u  £  Hk( fl)  n  /fo(fi)  using  the  space  Sp  C  H£{U)  of 
conforming  piecewise  polynomial  functions  of  degree  <  p  was  shown  to  satisfy 

IMI»‘(n)  <  C(€)p-{*-I)+c||u||tfnn)  (2.1) 

for  €  >  0  arbitrarily  small. 

The  proof  of  (2.1)  indicated  that  the  term  C(e)  could  grow  quickly  with  e  — ♦  0. 
Nevertheless,  computational  experience  suggested  that  (2.1)  holds  without  the  term 
e,  i.e. 


IMI/zqn)  <  Cj>  (lc  ‘’iMI/fUm  (2.2) 

Therefore,  the  first  problem  considered  by  us  was  to  solidify  the  foundations  of 
the  p-  and  h  —  p  versions  by  proving  optimal  estimates  like  (2.2)  for  the  expected 
rates  of  convergence.  These  results  are  summarized  below. 

(a)  In  [B2]  we  have  proven  optimal  approvimalion  results  using  the  p-version 
with  G,n  elements.  This  includes  the  proof  ol  (2.2)  for  the  case  when  the  solution 
is  smooth  and  a  separate  estimate  for  the  case  of  a  singular  solution. 


(b)  In  [B3j  similar  results  are  proven  for  the  k  —  p  version  using  quasiuniform 
meshes.  Numerical  experiments  are  described  which  show  the  higher  rates  of 
convergence  obtained  with  the  h  —  p  (and  p )  version  over  the  h-version. 

(c)  In  [Bli]  the  results  from  [B2]  are  extended  to  the  case  when  elements  pos¬ 
sessing  C”1-1  continuity  are  used.  These  are  necessary  for  problems  in  elasticity, 
like  those  involving  plates  and  shells.  Our  results  improve  upon  previous  work 
by  Katz  and  Wang  (for  m  —  2)  [A8j  and  by  Dorr  (for  general  m)  [A7]  in  two 
ways.  First,  like  (2.1),  the  results  in  [A8],  [A7]  are  optimal  up  to  an  arbitrary 
c  >  0,  while  ours  are  optimal.  Second,  in  order  to  prove  estimates  for  the  case 
of  elements  satisfying  a  homogeneous  essential  boundary  condition  in  H™( fl), 
the  authors  in  [A8]  and  [A7]  are  forced  to  make  an  assumption  involving  the 
interpolation  of  Sobolev  spaces  ([ A8]  eqn.  2.37  and  implicitly  in  [A7]  Theorem 
3.4).  We  are  able  to  prove  our  estimate  without  using  this  assumption  provided 
our  solution  lies  in  Hk{ ft),  k  >  2m  —  Since  we  have  a  separate  estimate  to 
treat  the  singular  portion  of  the  solution,  this  restriction  is  not  a  severe  one  in 
practice.  (Using  the  interpolation  assumption  allows  us  to  extend  our  results 
for  the  case  m  <  k  <  2m  —  j). 

2.2.  Approximation  of  Boundary  Conditions  by  the  h-  and  h-p  versions 

Previous  results  for  the  p-version  in  [A4]  and  [A7]  and  for  the  h  —  p  version  in  [A3]  did 
not  deal  with  the  case  of  inhomogeneous  essential  boundary  conditions.  Consequently, 
the  first  version  of  PROBE  assumed  that  essential  boundary  data  was  either  zero  or 
a  piecewise  polynomial  of  fixed  low  degree  on  the  grid  being  used.  As  a  result,  for 
general  inhomogeneous  conditions,  when  higher  degree  polynomials  were  used  inside, 
the  accuracy  of  the  approximation  on  the  boundary  did  not  increase  correspondingly, 
leading  to  a  loss  in  the  rate  of  convergence  in  this  case.  What  was  needed  was  a 
method  that  would  simultaneously  increase  the  degree  of  approximation  inside  fl  and 
on  the  boundary.  The  results  obtained  by  us  to  address  this  problem  are  summarized 
below. 

(a)  In  [B2j,  we  presented  and  analyzed  a  method  to  approximate  the  boundary 
data  based  on  a  projection  in  the  // 1  norm.  A  similar  method  for  the  h-p 
version  was  analyzed  in  [B3j.  Our  results  had  a  quick  practical  payofT,  since 
these  methods  were  incorporated  into  the  new  version  of  PROBE. 

(b)  hi  [B8|,  we  showed  that  the  If1  pro jer t  i<>M  in.-<\  not  work  ;is  well  for  problems 
with  rough  Dirichlet  data  on  the  boundary  An  alternate  method,  based  on  a 
projection  in  the  H *  norm  was  introduced  and  shown  to  be  more  robust.  This 
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method  generalizes  easily  to  systems  of  equations  and  can  be  implemented  using 
Fast  Fourier  Transforms. 

(c)  In  [B 13] ,  we  analyzed  a  family  of  methods  for  inhomogeneous  Dirichlet 
conditions  based  on  the  use  of  Jacobi  polynomials  with  varying  indices.  The 
H1  and  H 2  projections  analyzed  in  [B2]  and  [B8]  are  special  members  of  this 
family.  We  considered  the  comparison  of  various  projections  by  looking  at 
both  their  theoretical  and  computational  properties.  It  was  shown  numerically 
that  for  certain  problems,  the  use  of  an  incorrect  boundary  data  approximation 
method  could  lead  to  a  complete  lack  of  convergence  for  the  p-method. 

(d)  In  [B9],  we  considered  the  implementation  of  inhomogeneous  conditions 
for  elliptic  systems  of  equations,  particularly  “boundary  constraint  problems”. 
A  mixed-type  method  using  the  p-version  for  this  problem  is  described  and 
analyzed  in  the  above  reference. 

A  detailed  summary  of  the  results  in  Section  2.1  and  2.2  has  been  presented  by 
us  in  the  expository  paper  [B4] .  We  have  also  presented  some  of  these  results  in  [Bl] 
and  [B7j. 


2.3.  The  p-version  for  Mixed  Methods 


Mixed  finite  element  methods  are  useful  for  direct  approximation  of  physical  quanti¬ 
ties  of  interest  (which  may  not  be  the  primary  unknown  in  the  usual  formation  of  the 
problem).  They  have  obtained  an  enormous  amount  of  attention  in  the  literature  in 
the  context  of  the  classical  /i-version.  Unlike  standard  methods,  where  convergence 
only  depends  on  approximability  of  the  subspaces  used,  the  convergence  of  mixed 
methods  also  depends  upon  the  stability  of  the  subspaces.  Usually,  this  amounts  to 
ensuring  that  a  compatibility  condition  (the  Babuska-Brezzi  or  inf-sup  condition)  is 
satisfied  by  the  pair  of  surfaces  used. 

We  have  analyzed  the  p-version  formulations  of  three  mixed  methods  which  have 
been  previously  investigated  in  the  context  of  the  /i-version.  In  each  case,  the  ques¬ 
tion  of  approximability  is  settled  by  using  the  results  described  in  Section  2.1.  The 
question  of  stability  requires  a  different  approach. 


(a)  Our  first  result  concerns  the  solution  of  typical  elliptic  problems  like  the  Pois¬ 
son  equation  over  a  polygonal  domain  with  boundary  T  =  (J,n=i  T,.  We  may  write  this 
as  a  mixed  method  by  introducing  a  Lagrange  multiplier  4>  for  the  Dirichlet  boundary 
data.  It  may  be  shown  that  <t>  is  formally  oqoiv .* I >  m  !•»  i|i-  normal  'omponent  of  the 

flux  along  the  boundary  (i.e.  <j)  =  —  — ). 

dn 
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To  approximate  our  problem,  we  introduce  finite-dimensional  subspaces  Vp  C 
ff'(Sl)  ,  Sp  C  H  5(D,  parametrized  by  p,  the  degree  of  polynomials  used.  Here, 
functions  in  Vp  are  piecewise  polynomials  of  degree  <  p  on  a  fixed  grid  on  ft,  while 
functions  in  Sp  are  piecewise  polynomials  of  degree  <  p-  2  on  each  I\.  Then  we  seek 
a  pair  ( up,4>p )  £  Vp  x  SP  satisfying 


L 


Vup  •  Vvp  + 


/  vp<y 
Jr 

=  /  fv- 
Jr 

Vvp  £  Vp° 

[  upV 
r 

fe. 

Oi 

^ 

II 

V0P  £  Sp 

up(iV,) 

for  each  node  Ni  £  F 

(2.3) 

where  V^0  c  Vp  consists  of  those  functions  vanishing  at  each  TV,  £  I\ 

We  have  proven  that  if  u  £  Hk(Q)  with  k  >  then  the  error  in  u  satisfies  the 
optimal  estimate 


11“  ~  «Plltf*(n)  <  Cp  {k  1)||u||//*(n) 

(2.4) 

Moreover,  the  inf-sup  condition  takes  the  following  form: 

inf 

sup  /  v<j>  >  Cp  * 

J  r 

<t>esp 

v  £  Vp 

IWIfH(n  =  1 

IMU'fn)  =  1 

(2.5) 

which  leads  to  the  error  estimate 

■$(rj  <  Cp  (*+i)||</>||//*(r) 

(2.6) 

(Ilere  ||  •  ||„.(p,  :=  ||  •  ||„.(r.,) 

The  result  (2.4)  has  been  presented  by  us  in  [B9],  where  we  have  used  it  in  the 
analysis  of  the  boundary  constraint  problem.  The  results  (2.5)-(2.6)  have  not  been 
published  yet. 


(b)  An  alternative  mixed  formulation  of  the  Poisson  equation  involves  writing  it 
as  a  first  order  system  so  that  u  and  d  —  grad  u  are  the  unknowns  in  the  spaces  L2(ft) 
and  H  (div,H)  respectively.  The  most  well  known  finite  element  method  spaces  for 
this  method  are  the  Raviart-Thomas  spares  [.\"|  ! /,-((?).  //(,’  If  (div,  O) 
which  have  been  defined  for  arbitrary  polynomial  dogirc  p.  These  spares  have  been 
thoroughly  analyzed  in  the  context  of  the  /i-version.  In  particular,  it  has  been  shown 
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that  when  p  —  Po  is  fixed,  they  satisfy  the  inf-sup  condition  for  the  h  version  (see 
[A9]),  i.e. 

inf  sup  /rdiva«>c 

V  e  v£"  a  e  HZ' 

IMUJ(n)  =  1  il^H//(div,n)  =  ^ 
where  C  is  a  constant  independent  of  h  but  depending  upon  p0.  (2.7)  leads  to  optimal 
error  estimates  for  the  error  in  terms  of  h. 

Suppose  now  that  we  use  a  fixed  grid  of  mesh  size  h0  and  increase  p,  the  degree 
of  the  spaces  used  instead,  i.e.,  we  consider  the  pairs  V£(}  with  increasing  p. 

We  would  like  to  know  whether  the  method  based  on  this  p-type  extension  process  is 
stable  or  not. 

In  [B15],  we  have  analyzed  the  question  of  stability  for  arbitrary  combinations  of 
h  and  p.  In  particular,  we  have  shown  that  (2.7)  holds  for  the  p-version  as  well,  so 
that  C  is  a  constant  independent  of  p  and  h.  We  have  also  derived  estimates  that  are 
optimal  in  both  h  and  p  (up  to  an  arbitrary  e  >  0).  Our  results  indicate  in  addition 
that  the  use  of  the  /i-ve rsion  with  p  >  2  is  to  be  recommended  over  the  /i-version 
with  p  =  1. 

(c)  The  Brezzi-Douglas-Marini  elements  [A6]  are  an  alternative  to  the  Raviart- 
Thomas  elements.  We  have  shown  in  [B15j  that  the  inf-sup  condition  for  these  ele¬ 
ments  satisfies  (2.7)  where  for  arbitrary  e  >  0, 

C(p)  >  Cp-‘  (2.8) 

Here  C  depends  upon  (  but  is  independent  of  p  and  h.  This  once  again  leads  to 
optimal  estimates  (up  to  «  >  0)  in  both  p  and  h  [i.e.  for  the  h-,  p-,  and  h-p  versions] 

2.4.  The  p-  and  h-p  versions  for  Integral  Equations 

Many  problems  that  are  posed  originally  over  a  region  0  with  boundary  dfi  may  be 
re-formulated,  using  integral  equations,  as  problems  on  just  dU  alone.  We  believe  that 
the  p-  and  h-p  versions  will  prove  to  be  valuable  tools  in  solving  Galerkin  formulations 
of  integral  equations. 

(a)  In  [BIO),  we  have  analyzed  the  p-version  for  integral  equations  of  the  first  kind 
that  arise  form  the  Dirichlet  crack  problem  and  the  Neumann  crack  problem  in  two- 
dimensional  elasticity.  Our  results  also  apply  to  the  corresponding  screen  problems 
in  two-dimensional  acoustics.  The  solutions  fin  Mi--..  :,t,.  known  to  have 

an  r"  type  singularity,  which  leads  to  the  ml*-  *>l  <  on \<i  g<n<  <■  (in  the  energy 

nrm)  as  that  observed  using  the  /i-version.  We  have  also  analyzed  p-version  methods 
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for  some  other  integral  equations  arising  from  two  -  and  three-  dimensional  problems 
arising  in  elasticity  and  acoustics.  As  expected,  the  p- version  yields  higher  rates  of 
convergence  in  most  cases. 

(b)  In  [B 14],  we  extend  our  results  to  the  h-p  version.  We  consider  integral  equations 
on  polygonal  two-dimensional  domains,  for  which  the  solution  is  known  to  possess 
singularities  of  the  form  r°  where  a  is  general  and  depends  upon  the  geometry  of  the 
domain.  We  consider  the  case  with  mixed  boundary  conditions  and  derive  estimates 
that  take  into  consideration  both  the  degree  p  and  mesh-size  h  of  the  quasiuniform 
mesh  used. 

Our  results  provide  a  partial  theoretical  basis  to  some  computations  reported 
by  E.  Rank  on  the  p-version,  where  our  predicted  rate  of  convergence  was  observed 
experimentally  for  a  range  of  p  between  1  and  8.  (See  [B 14 ]  for  details). 

Summaries  of  the  above  results  may  be  found  in  [B6],  [B12], 

2.5.  An  implementation  of  the  p-  and  h-p  versions  for  MODULEF 

The  MODULEF  club,  created  by  INRIA,  Rocquencourt,  France  has  141  members 
worldwide  and  brings  together  universities  and  industrial  companies  with  the  goal 
of  designing  and  implementing  a  library  of  finite  element  modules.  The  p-  and  h-p 
versions  were  implemented  for  scalar  linear  elliptic  problems  during  an  extended  visit 
there,  in  the  form  of  a  module  called  ‘HP’.  This  module  uses  polynomials  of  degree 
1  <  P  <  8  combined  with  geometrically  refined  meshes  to  attain  approximations  with 
a  high  degree  of  accuracy.  A  feature  of  HP  allows  the  user  to  obtain  (at  low  additional 
cost)  a  sequence  of  solutions  for  various  p,  providing  an  effective  and  easy  method  to 
decide  whether  the  convergence  of  a  solution  is  acceptable  or  not.  (See  [B16]) 

HP  is  the  only  module  in  MODULEF  that  is  not  based  on  the  classical  h-version. 
The  availability  of  this  cede  to  the  large  body  of  MODULEF  users  is  expected  to 
greatly  enhance  the  awareness,  use  and  understanding  of  these  new  versions  and 
stimulate  further  research  in  this  field.  The  code  will  be  thoroughly  tested  before 
being  released  in  1989. 

2.6.  Numerical  methods  for  a  class  of  reaction-diffusion  problems 

We  consider  a  class  of  steady-state  semilinear  reaction-diffusion  problems  with  non- 
differentiablc  kinetics  whose  analytical  properties  have  received  considerable  attention 
in  the  literature,  see  e.g.  j  A5j.  More  precisely,  we  wish  to  solve  the  following  equations 
for  the  concentration  u: 

A u  -  A/(u(x))  x  e  D  (u  >  0  in  U)  (2.9) 
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u  =  1  xedn  (2.10) 

where  n  C  IRn  ( n  =  1,2,3)  is  a  bounded  domain  with  smooth  boundary  r?fi  and 
A  is  a  positive  constant  measuring  the  ratio  of  reaction  to  diffusion  rates.  We  are 
interested  in  the  case  of  pth  order  isothermal  reactions,  for  the  case  0  <  p  <  1,  where 
/  is  explicitly  given  by 


f{t)  =  tv  for  t  >  0 

=  0  for  t  <  0  (2.11) 

For  p  >  1  this  nonlinear  problem  may  be  approximated  by  several  methods. 
However,  when  p  <  1,  /  is  not  differentiable  at  the  origin  and  the  usual  techniques 
of  analysis  do  not  work. 

In  [B5],  we  have  considered  a  finite  element  approximation  for  (2.9)-(2.1l).  We 

have  shown  that  the  approximate  problem  has  a  unique  solution  vh,  that  when  linear 

functions  are  used,  vh  >  0,  and  that  the  following  error  estimate  can  be  obtained: 

ll«  -  v^llz/qn)  <  Ch  if  l-  <  p  <  1 

<  Ch2p  if  0  <  p  <  1  (2.12) 

2 

Numerical  experiments  reported  by  us  show,  however,  that  one  gets  0(/i)  convergence 
for  any  p  so  that  the  result  (2.12)  is  pessimistic  for  the  case  p  <  ]. 

We  also  analyze  a  finite  difference  scheme  for  (2.9)-(2.11),  proving  existence, 
uniqueness  and  convergence  of  the  approximation  using  the  theory  of  M-functions. 
Details  of  the  analysis  together  with  numerical  experiments  in  one-  and  two-dimensions 
may  be  found  in  [B5]. 
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